home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / b / b.lha / B / src / bint / b1nuR.c < prev    next >
C/C++ Source or Header  |  1988-11-24  |  6KB  |  285 lines
  1. /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
  2.  
  3. /*
  4.   $Header: b1nuR.c,v 1.4 85/08/22 16:51:49 timo Exp $
  5. */
  6.  
  7. /* Rational arithmetic */
  8.  
  9. #include "b.h"
  10. #include "b0con.h"
  11. #include "b1obj.h"
  12. #include "b1num.h"
  13. #include "b3err.h"
  14.  
  15. /* Length calculations used for fraction sizes: */
  16. #define Maxlen(u, v) \
  17.     (Roundsize(u) > Roundsize(v) ? Roundsize(u) : Roundsize(v))
  18. #define Sumlen(u, v) (Roundsize(u)+Roundsize(v))
  19. #define Difflen(u, v) (Roundsize(u)-Roundsize(v))
  20.  
  21. /* To shut off lint and other warnings: */
  22. #undef Copy
  23. #define Copy(x) ((integer)copy((value)(x)))
  24.  
  25. /* Globally used constants */
  26.  
  27. rational rat_zero;
  28. rational rat_half;
  29.  
  30. /* Make a normalized rational from two integers */
  31.  
  32. Visible rational mk_rat(x, y, len) integer x, y; int len; {
  33.     rational a;
  34.     integer u,v;
  35.  
  36.     if (y == int_0) {
  37.         if (interrupted)
  38.             return rat_zero;
  39.         syserr(MESS(1200, "mk_rat(x, y) with y=0"));
  40.     }
  41.  
  42.     if (x == int_0 && len <= 0) return (rational) Copy(rat_zero);
  43.  
  44.     if (Msd(y) < 0) {    /* interchange signs */
  45.         u = int_neg(x);
  46.         v = int_neg(y);
  47.     } else {
  48.         u = Copy(x);
  49.         v = Copy(y);
  50.     }
  51.  
  52.     a = (rational) grab_rat();
  53.     if (len > 0 && len+2 <= Maxintlet) Length(a) = -2 - len;
  54.  
  55.     if (u == int_0 || v == int_1) {
  56.         /* No simplification possible */
  57.         Numerator(a) = Copy(u);
  58.         Denominator(a) = int_1;
  59.     } else {
  60.         integer g, abs_u;
  61.  
  62.         if (Msd(u) < 0) abs_u = int_neg(u);
  63.         else abs_u = Copy(u);
  64.         g = int_gcd(abs_u, v);
  65.         release((value) abs_u);
  66.  
  67.         if (g != int_1) {
  68.             Numerator(a) = int_quot(u, g);
  69.             Denominator(a) = int_quot(v, g);
  70.         } else {
  71.             Numerator(a) = Copy(u);
  72.             Denominator(a) = Copy(v);
  73.         }
  74.         release((value) g);
  75.     }
  76.  
  77.     release((value) u); release((value) v);
  78.  
  79.     return a;
  80. }
  81.  
  82.  
  83. /* Arithmetic on rational numbers */
  84.  
  85. /* Shorthands: */
  86. #define N(u) Numerator(u)
  87. #define D(u) Denominator(u)
  88.  
  89. Visible rational rat_sum(u, v) register rational u, v; {
  90.     integer t1, t2, t3, t4;
  91.     rational a;
  92.  
  93.     t2= int_prod(N(u), D(v));
  94.     t3= int_prod(N(v), D(u));
  95.     t1= int_sum(t2, t3);
  96.     t4= int_prod(D(u), D(v));
  97.     a= mk_rat(t1, t4, Maxlen(u, v));
  98.     release((value) t1); release((value) t2);
  99.     release((value) t3); release((value) t4);
  100.  
  101.     return a;
  102. }
  103.  
  104.  
  105. Visible rational rat_diff(u, v) register rational u, v; {
  106.     integer t1, t2, t3, t4;
  107.     rational a;
  108.  
  109.     t2= int_prod(N(u), D(v));
  110.     t3= int_prod(N(v), D(u));
  111.     t1= int_diff(t2, t3);
  112.     t4= int_prod(D(u), D(v));
  113.     a= mk_rat(t1, t4, Maxlen(u, v));
  114.     release((value) t1); release((value) t2);
  115.     release((value) t3); release((value) t4);
  116.  
  117.     return a;
  118. }
  119.  
  120.  
  121. Visible rational rat_prod(u, v) register rational u, v; {
  122.     integer t1, t2;
  123.     rational a;
  124.  
  125.     t1= int_prod(N(u), N(v));
  126.     t2= int_prod(D(u), D(v));
  127.     a= mk_rat(t1, t2, Sumlen(u, v));
  128.     release((value) t1); release((value) t2);
  129.  
  130.     return a;
  131. }
  132.  
  133.  
  134. Visible rational rat_quot(u, v) register rational u, v; {
  135.     integer t1, t2;
  136.     rational a;
  137.  
  138.     if (Numerator(v) == int_0) {
  139.         error(MESS(1201, "in u/v, v is zero"));
  140.         return (rational) Copy(rat_zero);
  141.     }
  142.  
  143.     t1= int_prod(N(u), D(v));
  144.     t2= int_prod(D(u), N(v));
  145.     a= mk_rat(t1, t2, Difflen(u, v));
  146.     release((value) t1); release((value) t2);
  147.  
  148.     return a;
  149. }
  150.  
  151.  
  152. Visible rational rat_neg(u) register rational u; {
  153.     register rational a;
  154.  
  155.     /* Avoid a real subtraction from zero */
  156.  
  157.     if (Numerator(u) == int_0) return (rational) Copy(u);
  158.  
  159.     a = (rational) grab_rat();
  160.     N(a) = int_neg(N(u));
  161.     D(a) = Copy(D(u));
  162.     Length(a) = Length(u);
  163.  
  164.     return a;
  165. }
  166.  
  167.  
  168. /* Rational number to the integral power */
  169.  
  170. Visible rational rat_power(a, n) rational a; integer n; {
  171.     integer u, v, tu, tv, temp;
  172.  
  173.     if (n == int_0) return mk_rat(int_1, int_1, 0);
  174.  
  175.     if (Msd(n) < 0) {
  176.         if (Numerator(a) == int_0) {
  177.             error(MESS(1202, "in 0**n, n is negative"));
  178.             return (rational) Copy(a);
  179.         }
  180.         if (Msd(Numerator(a)) < 0) {
  181.             u= int_neg(Denominator(a));
  182.             v = int_neg(Numerator(a));
  183.         }
  184.         else {
  185.             u = Copy(Denominator(a));
  186.             v = Copy(Numerator(a));
  187.         }
  188.         n = int_neg(n);
  189.     } else {
  190.         if (Numerator(a) == int_0) return (rational) Copy(a);
  191.             /* To avoid necessary simplification later on */
  192.         u = Copy(Numerator(a));
  193.         v = Copy(Denominator(a));
  194.         n = Copy(n);
  195.     }
  196.  
  197.     tu = int_1;
  198.     tv = int_1;
  199.  
  200.     while (n != int_0 && !interrupted) {
  201.         if (Odd(Lsd(n))) {
  202.             if (u != int_1) {
  203.                 temp = tu;
  204.                 tu = int_prod(u, tu);
  205.                 release((value) temp);
  206.             }
  207.             if (v != int_1) {
  208.                 temp = tv;
  209.                 tv = int_prod(v, tv);
  210.                 release((value) temp);
  211.             }
  212.             if (n == int_1)
  213.                 break; /* Avoid useless last squaring */
  214.         }
  215.  
  216.         /* Square u, v */
  217.  
  218.         if (u != int_1) {
  219.             temp = u;
  220.             u = int_prod(u, u);
  221.             release((value) temp);
  222.         }
  223.         if (v != int_1) {
  224.             temp = v;
  225.             v = int_prod(v, v);
  226.             release((value) temp);
  227.         }
  228.  
  229.         n = int_half(n);
  230.     } /* while (n!=0) */
  231.  
  232.     release((value) n);
  233.     release((value) u);
  234.     release((value) v);
  235.     a = (rational) grab_rat();
  236.     Numerator(a) = tu;
  237.     Denominator(a) = tv;
  238.  
  239.     return a;
  240. }
  241.  
  242.  
  243. /* Compare two rational numbers */
  244.  
  245. Visible relation rat_comp(u, v) register rational u, v; {
  246.     int sd, su, sv;
  247.     integer nu, nv;
  248.  
  249.     /* 1. Compare pointers */
  250.     if (u == v || N(u) == N(v) && D(u) == D(v)) return 0;
  251.  
  252.     /* 2. Either zero? */
  253.     if (N(u) == int_0) return int_comp(int_0, N(v));
  254.     if (N(v) == int_0) return int_comp(N(u), int_0);
  255.  
  256.     /* 3. Compare signs */
  257.     su = Msd(N(u));
  258.     sv = Msd(N(v));
  259.     su = (su>0) - (su<0);
  260.     sv = (sv>0) - (sv<0);
  261.     if (su != sv) return su > sv ? 1 : -1;
  262.  
  263.     /* 4. Compute numerator of difference and return sign */
  264.     nu= int_prod(N(u), D(v));
  265.     nv= int_prod(N(v), D(u));
  266.     sd= int_comp(nu, nv);
  267.     release((value) nu); release((value) nv);
  268.     return sd;
  269. }
  270.  
  271. Visible Procedure rat_init() {
  272.     rat_zero = (rational) grab_rat();
  273.     Numerator(rat_zero) = int_0;
  274.     Denominator(rat_zero) = int_1;
  275.  
  276.     rat_half = (rational) grab_rat();
  277.     Numerator(rat_half) = int_1;
  278.     Denominator(rat_half) = int_2;
  279. }
  280.  
  281. Visible Procedure endrat() {
  282.     release((value) rat_zero);
  283.     release((value) rat_half);
  284. }
  285.